require(estimatr)

#### data preparation ####
FtF.data <- read.csv("FtF_data.csv")
online.data <- read.csv("online_data.csv")

# dichotomize social media use
FtF.data$social.media.d <- (FtF.data$social.media > 1) * 1
online.data$social.media.d <- (online.data$social.media > 1) * 1

# dichotomize popularity in neighborhood
FtF.data$neighborhood.d <- (FtF.data$neighborhood > 2) * 1
online.data$neighborhood.d <- (online.data$neighborhood > 2) * 1

## function for subgroup comparison of SDB
SDB.group.test <- function(y, x, data, treat, block, DQ) {
  temp.data.1 <- subset(data, data[, x] == 0)
  n.obs.1 <- nrow(temp.data.1)
  list.reg.1 <- lm_robust(formula(paste0(y, " ~ ", treat)), data = temp.data.1, 
                          fixed_effects = paste0("~ ", block))
  list.est.1 <- as.numeric(list.reg.1$coefficients)
  list.se.1 <- as.numeric(list.reg.1$std.error)
  DQ.est.1 <- mean(temp.data.1[, DQ], na.rm = TRUE)
  DQ.se.1 <- sqrt(DQ.est.1 * (1 - DQ.est.1) / sum(! is.na(temp.data.1[, DQ])))
  SDB.est.1 <- DQ.est.1 - list.est.1
  SDB.se.1 <- sqrt(list.se.1 ^ 2 + DQ.se.1 ^ 2)
  SDB.p.1 <- (1 - pnorm(abs(SDB.est.1 / SDB.se.1))) * 2
  temp.data.2 <- subset(data, data[, x] == 1)
  n.obs.2 <- nrow(temp.data.2)
  list.reg.2 <- lm_robust(formula(paste0(y, " ~ ", treat)), data = temp.data.2, 
                          fixed_effects = paste0("~ ", block))
  list.est.2 <- as.numeric(list.reg.2$coefficients)
  list.se.2 <- as.numeric(list.reg.2$std.error)
  DQ.est.2 <- mean(temp.data.2[, DQ], na.rm = TRUE)
  DQ.se.2 <- sqrt(DQ.est.2 * (1 - DQ.est.2) / sum(! is.na(temp.data.2[, DQ])))
  SDB.est.2 <- DQ.est.2 - list.est.2
  SDB.se.2 <- sqrt(list.se.2 ^ 2 + DQ.se.2 ^ 2)
  SDB.p.2 <- (1 - pnorm(abs(SDB.est.2 / SDB.se.2))) * 2
  SDB.diff.est <- SDB.est.2 - SDB.est.1
  SDB.diff.se <- sqrt(list.se.1 ^ 2 + DQ.se.1 ^ 2 + list.se.2 ^ 2 + DQ.se.2 ^ 2)
  SDB.diff.p <- (1 - pnorm(abs(SDB.diff.est / SDB.diff.se))) * 2
  out <- list(summary = c(SDB.est.1, SDB.est.2, SDB.diff.est, 
                          SDB.diff.est + qnorm(0.025) * SDB.diff.se, 
                          SDB.diff.est + qnorm(0.975) * SDB.diff.se, 
                          SDB.diff.p), 
              details = list(SDB.est.1 = SDB.est.1, SDB.se.1 = SDB.se.1, 
                             SDB.p.1 = SDB.p.1, n.obs.1 = n.obs.1, 
                             SDB.est.2 = SDB.est.2, SDB.se.2 = SDB.se.2, 
                             SDB.p.2 = SDB.p.2, n.obs.2 = n.obs.2, 
                             SDB.diff.est = SDB.diff.est, 
                             SDB.diff.se = SDB.diff.se, 
                             SDB.diff.p = SDB.diff.p))
  names(out$summary) <- c("SDB est (x = 0)", "SDB est (x = 1)", "SDB diff est", 
                          "SDB diff lower", "SDB diff upper", "SDB diff p")
  out
}

#### tests of Hypotheses 3-6 (FtF) ####
EJK.result.FtF <- SDB.group.test("list", "EJK", FtF.data, 
                                 "treatment", "block", "DQ")
round(EJK.result.FtF$summary, 3)

RT.result.FtF <- SDB.group.test("list", "RT", FtF.data, 
                                "treatment", "block", "DQ")
round(RT.result.FtF$summary, 3)

social.media.result.FtF <- SDB.group.test("list", "social.media.d", FtF.data, 
                                          "treatment", "block", "DQ")
round(social.media.result.FtF$summary, 3)

neighborhood.result.FtF <- SDB.group.test("list", "neighborhood.d", FtF.data, 
                                          "treatment", "block", "DQ")
round(neighborhood.result.FtF$summary, 3)

#### tests of Hypotheses 3-6 (online) ####
EJK.result.online <- SDB.group.test("list", "EJK", online.data, 
                                    "treatment", "block", "DQ")
round(EJK.result.online$summary, 3)

RT.result.online <- SDB.group.test("list", "RT", online.data, 
                                   "treatment", "block", "DQ")
round(RT.result.online$summary, 3)

social.media.result.online <- SDB.group.test("list", "social.media.d", online.data, 
                                             "treatment", "block", "DQ")
round(social.media.result.online$summary, 3)

neighborhood.result.online <- SDB.group.test("list", "neighborhood.d", online.data, 
                                             "treatment", "block", "DQ")
round(neighborhood.result.online$summary, 3)

#### Figure 3 ####
cairo_pdf("Figure_3.pdf", width = 5, height = 7, pointsize = 9)
layout(matrix(1:2, 2, 1))
par(mar = c(4, 0, 3, 1), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-2.5, 1), ylim = c(0, 16), 
     main = "Face-to-face survey", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(-0.8, c(1:3, 5:7, 9:11, 13:15), 1.05, c(1:3, 5:7, 9:11, 13:15), 
         lty = 3, col = "gray")
abline(v = seq(-0.75, 1, 0.25), col = "gray")
points(EJK.result.FtF$details$SDB.est.1, 15, pch = 4)
segments(EJK.result.FtF$details$SDB.est.1 + 
           qnorm(0.025) * EJK.result.FtF$details$SDB.se.1, 15, 
         EJK.result.FtF$details$SDB.est.1 + 
           qnorm(0.975) * EJK.result.FtF$details$SDB.se.1, 15)
points(EJK.result.FtF$details$SDB.est.2, 14, pch = 4)
segments(EJK.result.FtF$details$SDB.est.2 + 
           qnorm(0.025) * EJK.result.FtF$details$SDB.se.2, 14, 
         EJK.result.FtF$details$SDB.est.2 + 
           qnorm(0.975) * EJK.result.FtF$details$SDB.se.2, 14)
points(EJK.result.FtF$details$SDB.diff.est, 13, pch = 19)
segments(EJK.result.FtF$details$SDB.diff.est + 
           qnorm(0.025) * EJK.result.FtF$details$SDB.diff.se, 13, 
         EJK.result.FtF$details$SDB.diff.est + 
           qnorm(0.975) * EJK.result.FtF$details$SDB.diff.se, 13)
text(EJK.result.FtF$details$SDB.diff.est, 13, 
     sprintf("%5.3f", EJK.result.FtF$details$SDB.diff.est), 
     cex = 0.8, pos = 1)
points(RT.result.FtF$details$SDB.est.1, 11, pch = 4)
segments(RT.result.FtF$details$SDB.est.1 + 
           qnorm(0.025) * RT.result.FtF$details$SDB.se.1, 11, 
         RT.result.FtF$details$SDB.est.1 + 
           qnorm(0.975) * RT.result.FtF$details$SDB.se.1, 11)
points(RT.result.FtF$details$SDB.est.2, 10, pch = 4)
segments(RT.result.FtF$details$SDB.est.2 + 
           qnorm(0.025) * RT.result.FtF$details$SDB.se.2, 10, 
         RT.result.FtF$details$SDB.est.2 + 
           qnorm(0.975) * RT.result.FtF$details$SDB.se.2, 10)
points(RT.result.FtF$details$SDB.diff.est, 9, pch = 19)
segments(RT.result.FtF$details$SDB.diff.est + 
           qnorm(0.025) * RT.result.FtF$details$SDB.diff.se, 9, 
         RT.result.FtF$details$SDB.diff.est + 
           qnorm(0.975) * RT.result.FtF$details$SDB.diff.se, 9)
text(RT.result.FtF$details$SDB.diff.est, 9, 
     sprintf("%5.3f", RT.result.FtF$details$SDB.diff.est), 
     cex = 0.8, pos = 1)
points(social.media.result.FtF$details$SDB.est.1, 7, pch = 4)
segments(social.media.result.FtF$details$SDB.est.1 + 
           qnorm(0.025) * social.media.result.FtF$details$SDB.se.1, 7, 
         social.media.result.FtF$details$SDB.est.1 + 
           qnorm(0.975) * social.media.result.FtF$details$SDB.se.1, 7)
points(social.media.result.FtF$details$SDB.est.2, 6, pch = 4)
segments(social.media.result.FtF$details$SDB.est.2 + 
           qnorm(0.025) * social.media.result.FtF$details$SDB.se.2, 6, 
         social.media.result.FtF$details$SDB.est.2 + 
           qnorm(0.975) * social.media.result.FtF$details$SDB.se.2, 6)
points(social.media.result.FtF$details$SDB.diff.est, 5, pch = 19)
segments(social.media.result.FtF$details$SDB.diff.est + 
           qnorm(0.025) * social.media.result.FtF$details$SDB.diff.se, 5, 
         social.media.result.FtF$details$SDB.diff.est + 
           qnorm(0.975) * social.media.result.FtF$details$SDB.diff.se, 5)
text(social.media.result.FtF$details$SDB.diff.est, 5, 
     sprintf("%5.3f", social.media.result.FtF$details$SDB.diff.est), 
     cex = 0.8, pos = 1)
points(neighborhood.result.FtF$details$SDB.est.1, 3, pch = 4)
segments(neighborhood.result.FtF$details$SDB.est.1 + 
           qnorm(0.025) * neighborhood.result.FtF$details$SDB.se.1, 3, 
         neighborhood.result.FtF$details$SDB.est.1 + 
           qnorm(0.975) * neighborhood.result.FtF$details$SDB.se.1, 3)
points(neighborhood.result.FtF$details$SDB.est.2, 2, pch = 4)
segments(neighborhood.result.FtF$details$SDB.est.2 + 
           qnorm(0.025) * neighborhood.result.FtF$details$SDB.se.2, 2, 
         neighborhood.result.FtF$details$SDB.est.2 + 
           qnorm(0.975) * neighborhood.result.FtF$details$SDB.se.2, 2)
points(neighborhood.result.FtF$details$SDB.diff.est, 1, pch = 19)
segments(neighborhood.result.FtF$details$SDB.diff.est + 
           qnorm(0.025) * neighborhood.result.FtF$details$SDB.diff.se, 1, 
         neighborhood.result.FtF$details$SDB.diff.est + 
           qnorm(0.975) * neighborhood.result.FtF$details$SDB.diff.se, 1)
text(neighborhood.result.FtF$details$SDB.diff.est, 1, 
     sprintf("%5.3f", neighborhood.result.FtF$details$SDB.diff.est), 
     cex = 0.8, pos = 1)
text(-2.5, 16, "Extrajudicial killing", pos = 4, font = 2)
text(-2.35, 15, paste0("Do not know a victim (N=", 
                       format(EJK.result.FtF$details$n.obs.1, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 14, paste0("Know a victim (N=", 
                       format(EJK.result.FtF$details$n.obs.2, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 13, "Difference b/w groups", pos = 4)
text(-2.5, 12, "Red-tagging", pos = 4, font = 2)
text(-2.35, 11, paste0("Do not know a victim (N=", 
                       format(RT.result.FtF$details$n.obs.1, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 10, paste0("Know a victim (N=", 
                       format(RT.result.FtF$details$n.obs.2, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 9, "Difference b/w groups", pos = 4)
text(-2.5, 8, "Social media use", pos = 4, font = 2)
text(-2.35, 7, paste0("60 minutes or less (N=", 
                      format(social.media.result.FtF$details$n.obs.1, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 6, paste0("More than 60 minutes (N=", 
                      format(social.media.result.FtF$details$n.obs.2, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 5, "Difference b/w groups", pos = 4)
text(-2.5, 4, "Popularity in neighborhood", pos = 4, font = 2)
text(-2.35, 3, paste0("Dissatisfied (N=", 
                      format(neighborhood.result.FtF$details$n.obs.1, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 2, paste0("Satisfied (N=", 
                      format(neighborhood.result.FtF$details$n.obs.2, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 1, "Difference b/w groups", pos = 4)
axis(1, seq(-0.75, 1, 0.25), lwd = 0.5, labels = NA)
mtext(sprintf("%4.2f", seq(-0.75, 1, 0.25)), 
      side = 1, line = 1, at = seq(-0.75, 1, 0.25))
mtext("Magnitude of SDB (or its difference)", side = 1, line = 2.5, at = 0.125)
plot(NULL, NULL, bty = "n", xlim = c(-2.5, 1), ylim = c(0, 16), 
     main = "Online survey", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(-0.8, c(1:3, 5:7, 9:11, 13:15), 1.05, c(1:3, 5:7, 9:11, 13:15), 
         lty = 3, col = "gray")
abline(v = seq(-0.75, 1, 0.25), col = "gray")
points(EJK.result.online$details$SDB.est.1, 15, pch = 4)
segments(EJK.result.online$details$SDB.est.1 + 
           qnorm(0.025) * EJK.result.online$details$SDB.se.1, 15, 
         EJK.result.online$details$SDB.est.1 + 
           qnorm(0.975) * EJK.result.online$details$SDB.se.1, 15)
points(EJK.result.online$details$SDB.est.2, 14, pch = 4)
segments(EJK.result.online$details$SDB.est.2 + 
           qnorm(0.025) * EJK.result.online$details$SDB.se.2, 14, 
         EJK.result.online$details$SDB.est.2 + 
           qnorm(0.975) * EJK.result.online$details$SDB.se.2, 14)
points(EJK.result.online$details$SDB.diff.est, 13, pch = 19)
segments(EJK.result.online$details$SDB.diff.est + 
           qnorm(0.025) * EJK.result.online$details$SDB.diff.se, 13, 
         EJK.result.online$details$SDB.diff.est + 
           qnorm(0.975) * EJK.result.online$details$SDB.diff.se, 13)
text(EJK.result.online$details$SDB.diff.est, 13, 
     sprintf("%5.3f", EJK.result.online$details$SDB.diff.est), 
     cex = 0.8, pos = 1)
points(RT.result.online$details$SDB.est.1, 11, pch = 4)
segments(RT.result.online$details$SDB.est.1 + 
           qnorm(0.025) * RT.result.online$details$SDB.se.1, 11, 
         RT.result.online$details$SDB.est.1 + 
           qnorm(0.975) * RT.result.online$details$SDB.se.1, 11)
points(RT.result.online$details$SDB.est.2, 10, pch = 4)
segments(RT.result.online$details$SDB.est.2 + 
           qnorm(0.025) * RT.result.online$details$SDB.se.2, 10, 
         RT.result.online$details$SDB.est.2 + 
           qnorm(0.975) * RT.result.online$details$SDB.se.2, 10)
points(RT.result.online$details$SDB.diff.est, 9, pch = 19)
segments(RT.result.online$details$SDB.diff.est + 
           qnorm(0.025) * RT.result.online$details$SDB.diff.se, 9, 
         RT.result.online$details$SDB.diff.est + 
           qnorm(0.975) * RT.result.online$details$SDB.diff.se, 9)
text(RT.result.online$details$SDB.diff.est, 9, 
     sprintf("%5.3f", RT.result.online$details$SDB.diff.est), 
     cex = 0.8, pos = 1)
points(social.media.result.online$details$SDB.est.1, 7, pch = 4)
segments(social.media.result.online$details$SDB.est.1 + 
           qnorm(0.025) * social.media.result.online$details$SDB.se.1, 7, 
         social.media.result.online$details$SDB.est.1 + 
           qnorm(0.975) * social.media.result.online$details$SDB.se.1, 7)
points(social.media.result.online$details$SDB.est.2, 6, pch = 4)
segments(social.media.result.online$details$SDB.est.2 + 
           qnorm(0.025) * social.media.result.online$details$SDB.se.2, 6, 
         social.media.result.online$details$SDB.est.2 + 
           qnorm(0.975) * social.media.result.online$details$SDB.se.2, 6)
points(social.media.result.online$details$SDB.diff.est, 5, pch = 19)
segments(social.media.result.online$details$SDB.diff.est + 
           qnorm(0.025) * social.media.result.online$details$SDB.diff.se, 5, 
         social.media.result.online$details$SDB.diff.est + 
           qnorm(0.975) * social.media.result.online$details$SDB.diff.se, 5)
text(social.media.result.online$details$SDB.diff.est, 5, 
     sprintf("%5.3f", social.media.result.online$details$SDB.diff.est), 
     cex = 0.8, pos = 1)
points(neighborhood.result.online$details$SDB.est.1, 3, pch = 4)
segments(neighborhood.result.online$details$SDB.est.1 + 
           qnorm(0.025) * neighborhood.result.online$details$SDB.se.1, 3, 
         neighborhood.result.online$details$SDB.est.1 + 
           qnorm(0.975) * neighborhood.result.online$details$SDB.se.1, 3)
points(neighborhood.result.online$details$SDB.est.2, 2, pch = 4)
segments(neighborhood.result.online$details$SDB.est.2 + 
           qnorm(0.025) * neighborhood.result.online$details$SDB.se.2, 2, 
         neighborhood.result.online$details$SDB.est.2 + 
           qnorm(0.975) * neighborhood.result.online$details$SDB.se.2, 2)
points(neighborhood.result.online$details$SDB.diff.est, 1, pch = 19)
segments(neighborhood.result.online$details$SDB.diff.est + 
           qnorm(0.025) * neighborhood.result.online$details$SDB.diff.se, 1, 
         neighborhood.result.online$details$SDB.diff.est + 
           qnorm(0.975) * neighborhood.result.online$details$SDB.diff.se, 1)
text(neighborhood.result.online$details$SDB.diff.est, 1, 
     sprintf("%5.3f", neighborhood.result.online$details$SDB.diff.est), 
     cex = 0.8, pos = 1)
text(-2.5, 16, "Extrajudicial killing", pos = 4, font = 2)
text(-2.35, 15, paste0("Do not know a victim (N=", 
                       format(EJK.result.online$details$n.obs.1, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 14, paste0("Know a victim (N=", 
                       format(EJK.result.online$details$n.obs.2, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 13, "Difference b/w groups", pos = 4)
text(-2.5, 12, "Red-tagging", pos = 4, font = 2)
text(-2.35, 11, paste0("Do not know a victim (N=", 
                       format(RT.result.online$details$n.obs.1, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 10, paste0("Know a victim (N=", 
                       format(RT.result.online$details$n.obs.2, big.mark = ","), 
                       ")"), pos = 4)
text(-2.35, 9, "Difference b/w groups", pos = 4)
text(-2.5, 8, "Social media use", pos = 4, font = 2)
text(-2.35, 7, paste0("60 minutes or less (N=", 
                      format(social.media.result.online$details$n.obs.1, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 6, paste0("More than 60 minutes (N=", 
                      format(social.media.result.online$details$n.obs.2, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 5, "Difference b/w groups", pos = 4)
text(-2.5, 4, "Popularity in neighborhood", pos = 4, font = 2)
text(-2.35, 3, paste0("Dissatisfied (N=", 
                      format(neighborhood.result.online$details$n.obs.1, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 2, paste0("Satisfied (N=", 
                      format(neighborhood.result.online$details$n.obs.2, big.mark = ","), 
                      ")"), pos = 4)
text(-2.35, 1, "Difference b/w groups", pos = 4)
axis(1, seq(-0.75, 1, 0.25), lwd = 0.5, labels = NA)
mtext(sprintf("%4.2f", seq(-0.75, 1, 0.25)), 
      side = 1, line = 1, at = seq(-0.75, 1, 0.25))
mtext("Magnitude of SDB (or its difference)", side = 1, line = 2.5, at = 0.125)
dev.off()